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A space-marching finite-volume method employing a nonorthogona! coordinate system and using a split dif- 
ferencing scheme for calculating steady supersonic flow over aerodynamic shapes is presented. It is a second - 
order-accurate mixed explicit-implicit procedure that solves the inviscid adiabatic and nondiffusive equations for 
chemically reacting flow in integral conservation-law form. The relationship between the finite-volume and dif- 
ferential forms of the equations Is examined and the relative merits of each discussed. The method admits initial 
Cauchy data situated on any arbitrary^ surface and integrates them forward along a general curvilinear coor- 
dinate, distorting and deforming the surface as it advances. The chemical kinetics term is split from the con- 
vective terms which are themselves dimensionally split, thereby freeing the fluid operators from the restricted 
step size imposed by the chemical reactions and increasing the computational efficiency. The accuracy of this 
splitting technique is analyzed, a sufficient stability criterion is established, a representative flow computation is 
discussed, and some comparisons are made with another method. 


Introduction 

C ONTINUING interest in and the usefulness of more 
sophisticated re-entry vehicles than those currently 
developed has culminated in the decision to design and con- 
struct a reusable space shuttle orbiter. This decision has rekin- 
dled the need for computing flows with finite-rate chemical 
reactions. The calculation of such flows is extremely useful to 
the vehicle designer as a source of information relative to the 
prediction of heat-transfer rates, boundary-layer effects, and 
aerodynamic loads acting on the aircraft. At moderate super- 
sonic velocities and lower altitudes the airflow remains in 
chemical equilibrium, and design information can be obtained 
from the flow calculations as well as wind-tunnel testing of 
appropriately scaled models. At higher velocities and 
altitudes, however, significant nonequilibrium effects are ex- 
pected to persist over a large area of the vehicle; in addition to 
affecting the surface heating rate, the dissociating flow 
produces atomic oxygen which might enhance the oxidation 
and erosion process at the surface. Complex chemically 
reacting flow phenomena cannot be scaled, and since current 
test facilities are incapable of carrying out full-scale ex- 
periments in this regime the designer can only turn to the com- 
puter solution for a realistic description of the flowfield at ac- 
tual flight conditions. 

For some time numerous procedures’"* for the calculation 
of one- and two-dimensional inviscid nonequilibrium flow 
have been known. Our discussion, however, is restricted to 
those techniques that take advantage of the hyperbolic charac- 
ter of the steady governing equations to integrate initial 
Cauchy data in a coordinate direction along which the local 
velocity component is everywhere supersonic. The method of 
characteristics has been applied recently to reacting flow in 
three-dimensional nozzles^ as well as around general-shaped 
bodies^; more recently a method for calculating reacting flow 
around an arbitrary body using a finite-difference procedure 
was introduced by Davy and Reinhardt. ’’ 

Common to all these existing techniques is the requirement 
that the initial data lie on a surface symmetrically oriented 
with respect to the marching coordinate, and the method of 
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Davy and Reinhardt even further stipulates that this initial 
surface be a plane normal to the body axis and that it advance 
downstream and undistorted. This condition simplifies the 
coordinate geometry and difference equations but severely 
restricts the choice of acceptable initial data and ignores the 
orientation of their cones of influence. For flow of a perfect 
gas Rizzi, et al.^ have devised a more general procedure that 
admits initial conditions situated on any arbitrary surface and 
integrates them forward along a general curvilinear coor- 
dinate that nearly conforms to the local streamline, rotating 
and distorting the initial data surface as it advances. Adopting 
that generalized marching method, this paper presents a 
method for computing the steady flow of an air mixture of 
dissociating and chemically nonequilibrated gases past a 
three-dimensional body. The theoretical description of the 
flow is idealized to the extent that diffusion, heat conduction, 
and viscous effects are considered negligible. 

In addition to increasing the number of equations needed to 
describe it, the chemical phenomena further complicate the 
numerical solution by introducing a second intrinsic rate 
process requiring a step size that can be quite different from 
that for convection alone. This paper proposes an efficient 
way of handling these two inherently different rates by 
developing finite-difference operators that split one spatial 
dimension from another as well as from the chemical kinetics 
term and that solve the governing equations to second-order 
accuracy. The convective operators use an explicit scheme, 
whereas the chemistry operator is implicit. An analysis of the 
accuracy of these split operators is carried out, and the com- 
putational efficiency gained by them is discussed. 

The initial data surface advancing downstream determines 
the mesh network for the next step, which is aligned with both 
the bow shock wave and the body surface. This alignment 
simplifies imposition of the boundary condition at both the 
shock and body, and ensures consistency with the interior 
flow at these points. 

Continuum Model 

The description of high-temperature airflow past a vehicle 
involves the calculation not only of the bulk-flow properties 
of every fluid element, but also the chemical phenomena that 
the individual gas species comprising the mixture are sub- 
jected to within the element. The chemical kinetics of air in 
the temperature range considered here, mainly entails five 
species -O2, N2, N, O, and NO -entering in the three 
dissociation reactions and the two bimolecular exchange reac- 
tions written in Table 1 . In these reactions M is the third-body 
collision partner that can be any of the five species in the mix- 
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Table 1 Chemical reactions 


f 

Chemical equation 

I 

02+M-20-hM 

2 

N2+M-2N-t-M 

3 

NO + M^^N + O + M 

4 

O 2 +N-NO + O 

5 

N2+0^N0-^N 


ture. The equations of motion for reacting flow must ef- 
fectively combine the microscopic formulation of chemical 
kinetics with the macroscopic equations of gas dynamics 
governing the bulk-flow properties. The appropriate descrip- 
tion for inviscid and nondiffusive flow of a mixture of react- 
ing gases in chemical nonequilibrium is the system of con- 
tinuum equations that incorporate these two diverse aspects 
into an overall macroscopic formulation. They are presented 
here without further discussion. Readers interested in their 
development and justification are referred to the textbook by 
Vincenti and Krueger. ^ 

The continuity equation for the mixture of density p is 

divpv = 0 (la) 

and the momentum equation takes the form 

div (pw-Hp /) (lb) 

where p signifies the macroscopic pressure, v the bulk 
velocity, and / is the identity tensor. The energy equation with 
enthalpy h may be written 

diwlpvih+Viv^)] =0 (1c) 

and can be integrated at once to obtain /?+ = constant 
along a streamline. For flow with no spatial variation in the 
freestream, the only case treated here, the constant is the same 
for all streamlines and this relation becomes 

h+Viv2=h, (Id) 

and hi equals the total enthalpy in the freestream and is a 
global constant. The macroscopic chemical kinetics equation 
for each of the five species is 

divpCfV = o:,>(p,r,C/,..,,Cj) (le) 

where T denotes the bulk temperature of the mixture, Cf 

represents the mass fraction p(/p, and stands for the 
chemical source function of species f, i.e., the mass rate of 
production of species f per unit volume by chemical reactions. 
Clearly, by definition, the condition 

f=i 

must hold. Vincenti and Krueger^ develop a suitable ex- 
pression for from principles of chemical kinetics, and the 
required constants and coefficients that we use are those given 
by Refs. 7 and 9. The system is completed by an equation of 
state in two forms, one relating enthalpy h = (p, p, ) and the 
other temperature T=T(p, k, C() to density, pressure, and 
species concentration, the particular forms of which are con- 
sidered next. 

In the uhdisTurbed flow, air is considered to be a mixture of 
23.3^0 oxygen and 76.7% nitrogen by mass, and in the cases 
dealt with here its temperature behind the bow shock varies 
typically between 2000 and 5000 K, seldom exceeding 6000 K. 
The mixture of pure air is assumed to remain trans- 
lationally and rotationally fully equilibrated while the 
vibrational excitation of nitric oxide and molecular oxygen 



Fig. 1 Generalized coordinates and flowfield discretization for 
space-marching finite-volume method. 

and nitrogen is approximated by the **half-excited” model of 
Lighthill, Rather than introducing nonequilibrium 
vibrational excitation, the vibrational mode is assumed to be 
excited to the extent that its energy is always half of the fully 
excited equilibrium value attained at high temperatures. The 
validity of the Lighthill model for the present application is 
given in Ref. 6. 

The five component species are thermally perfect gases in 
thermal equilibrium. With the above assumptions the state 
equation may be written explicitly as 

(2a) 

y-l P 

where h? represents the net heat per unit mass evolved during 
the formation of species f and 

r=p/p(Rj] (c./are,,) (2b) 

f=/ 

where the ratio of specific heats y is 

— (^o ) ^ 

^(^02 ^/v) 

and is the molecular weight of species f and (R is the 
universal gas constant. At this point note that the fluid and 
thermodynamic state of the gaseous system is determined 
completely from the calculation of p, q, v, and p satisfying 
the system of Eqs. (1) and (2) under specific boundary and 
initial conditions for a given problem. 

Differential and Integral Formulation 

For the method to admit initial Cauchy data situated on an 
arbitrary surface and march them in roughly the streamline 
direction, generalized coordinates must be used. The nonor- 
thogonal curvilinear system with contravariant coordinates 
and covariant base vectors gi (i= 1,2,3) is introduced 
where and lie in the initial data surface and x^ is in the 
approximate direction of the streamlinesj (Fig. 1). A 
reciprocal set of directions is associated with this system and 
is called contravariant field vectors grad x"' (m = 1,2,3) 
defined as the gradients of the curvilinear coordinates and 
related to the base vectors by =5 -"(Kronecker delta). 
Between these coordinates and a rectangular Cartesian system 
Zni with unit base vectors exists the functional relationship 

x^=x^(Zi,Z2,Z3) 1,2,3 (3) 

JTo minimize confusion, exponents are not used on vectors, coor- 
dinates, or vector components; superscripts denote only the con- 
travariem direction of the vectors and their components. 
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The governing Eqs. (I) div 3 = div {Fg] -\-Hg 2 + Ugs) =Q 
written in divergence form with respect to this curvilinear 
system is 


d(g'^U) , d(g'^H) , d{g''-F) _ ^ 
dx^ dx^ dx> 

where 

r 1 


(4) 


0) 






0 



Fig. 2 Formation of typical 
volume cell by the surfaces 
= constant. 


U= 


pu^Wi ~^p{dx^/dZj) 


fi(U) = 


0 


PU^W2 +p(dx^ /dZ2) 


0 


pu^Wj p(dx^ /dzj) 


0 


can be determined numerically at each mesh cell quickly by 
simple algebraic expressions: therein lies its advantage. 

Although some computational fluid dynamicists have used 
integral forms in solving flow problems, we are aware of no 


pCiU^ 


pCfU ' 

pu^ 


pu^ 

pu^Wi-^p(dx^/dZj ) 

F{U) = 

pu^ Wf +p(dx^ /dZ2) 

pU^W2-i-p(dx^ /dZ2) j 


pu^ W2-^p{dx^ /dZ2) 

pU^Wj+p(dx^/dZj) 


_pu^Wj-\-p(dx'/dZj) _ 


and g'" is the Jacobian d(zj, Zj»)/d{x\ x\ x, ) and the 
flow velocity is 

j 3 

i-I m=I 

A more physical interpretation of the contravarient coor- 
dinates and covariant directions can be achieved by using 
Gauss’ theorem to cast Eq. (4) into integral or, as we call it, 
finite-volume form 

+ |Fd5=|j Jfid vol (5) 

where 5^=5^'" is one of the six surfaces (with normal 
g^^) of the hexahedron formed by the coordinate surfaces 
= constant as illustrated in Fig. 2. Equations (4) and (5) are 
mathematically equivalent, but each has its particular ad- 
vantage. The numerical analysis of a method’s accuracy and 
stability is easier to execute with the partial differential form 
(4), whereas the actual computations are more conveniently 
carried out in the form (5); the latter is true because use of 
Eqs. (4) requires specifying the coordinate transformation (3) 
and calculating either analytically or numerically the Jacobian 
and the term dx’/3z,„ at each cell, which for practical 
body shapes may be time consuming and awkward. These two 
quantities, however, have a geometrical interpretation that 
becomes obvious in the finite-volume form. For example, 
dx^/dZm=S‘ simply the direction cosine between con- 

travarient direction / and Cartesian direction m and may be 
calculated from simple vector relations without having to per- 
form any partial differentiation. Similarly, the areas and 
volume which appear in Eq. (5), and are evaluated by dif- 
ferential expressions such as ds^ ^g '"^ dx^ dx^ and d vol = 
g- dx^ dx^ dx^, may also be calculated from geometrical 
principles. 1 The integral or finite-volume form in effect 
replaces explicit use of these functions with the products of 
vector quantities that have a definite physical significance and 

§Note that expressing the vector ^ = Fg/ Hgj Ugj involving 
products of mixed velocity components u' and of both bases 
g/ and a,„ leaves Eq. (4) in strong conservation form. See 8 and 11 
for further discussion of this point. 

^The area of each quadrilateral face of the hexahedron is one half 
of the vector cross product of the two diagonals on that face. The 
hexahedron itself is composed of five tetrahedra whose volumes are 
computed by a simple algebraic formula. 


one else who has shown the equivalence of the integral and 
curvilinear tensor differential forms, probably because of the 
classical bias for orthogonal coordinates. However, we em- 
phasize that these two formulations are identical, as of course 
they must be since they have the same governing equations; 
consequently, the accuracy and stability criteria of a chosen 
difference scheme when applied to one are exactly the same 
criteria when applied to the other. The only difference be- 
tween the two forms is the way in which the coordinate 
transformation terms are evaluated. The differential form is 
very useful for establishing the physical and numerical 
domains of dependence which determine the condition for 
numerical stability. It is also vital to the analysis of the ac- 
curacy of a given numerical procedure, but it is more con- 
venient to carry out the actual computations using the 
geometrical interpretation of the finite-volume form. 

Numerical Procedure 

Finite-difference approximations to the governing 
equations are used to advance the solution in thex^ direction 
from specified initial data. Split difference operators are con- 
structed in this section to solve the finite-volume formulation 
in Eq. (5) for computational cells like the one depicted in Fig. 
2 . 

Difference Operators 

The essential concept of the splitting technique is to “split” 
the governing equation into a series of simpler component 
equations, so that when solved sequentially in intermediate 
fractional steps they approximate the original equation to 
some order of accuracy. For nonreacting flow, MacCormack 
and Paullay,'" MacCormack and Warming, and Rizzi, et 
al.,^ among many others, have shown how this concept leads 
to splitting one spatial dimension from another; Rizzi and 
Bailey have extended the concept to the tirne-dependent 
equations of chemically reacting flow, and, in addition, have 
split the chemical production term from the spatial ones 
which are themselves dimensionally split. We apply this idea 
to the steady governing Eq. (5) and break them into three dif- 
ference operators 


^ -FU^,S, ) 


UjT = 


J,k + J^k + i ~ f^j,k ^k) 


(6a) 
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(6c) 

where \ = 2k=V3 are the fractional step numbers, and the af- 
fixal notation now no longer refers to the contravarient and 
covariam directions but rather to the geometrical location in 
the difference mesh. The subscripts j and k denote the discrete 
location of the cell j,k along the axes and x\ while 
Sj(Ax ^ ) and 5^^ {Ax - ) are the areas of its sides facing these 
directions. The superscripts denote the position 


= ^ Ax^n) 

n 


to which the solution has been advanced during each cycle of 
the operators, and S" is the area facing x^. The mean values 
Uj f^ of the flow variables in the cell at step n are defined by 


where 


ith 


o \dn/ 


dp / f>,c^ 


is the square of the local frozen speed of sound, y* = V'g’/ 
(g") is the component of v along the g‘ direc- 

tion, and 


g= — 2j/ V-^ cos o) 

is the component of v lying in the plane defined by and g^ 
and cos a=g^ *,g^/(g ”g”) and 

j 

g" = 2^ (dx'/dz„, ) (dx'/dz ,„ ) 

m = J 

A similar relation determining 6^ also exists for the L 2 
operator. The chemical rate operator L^. is implicit and a 
theoretical analysis using a linearized production term for 
A assures that it is unconditionally stable and hence may 
be controlled by either of the two explicit step sizes. The 
sequence of operators 

( 10 ) 


and —F( ) and = //( Ujj^ ) is the usual condensed 

notation. The matrix Qj of the fifth order is composed of 
elements )] while I is the identity 

matrix, and lastly vol^^^. (Ajc^)) is the volume of the 
hexahedron enclosed by the surfaces + + and 

^ This set of equations can be written more conveniently 
in operator notation. Let (Ax^) denote the operation per- 
formed by the set of Eqs. (6a) in advancing the solution from 
t/;* to UlV=L,(Ax^’)Ul,) and let L, (Ax’) 

and L,. (Ax'^) be similarly defined by Eqs. (6b) and (6c). 
Analogous to the familiar MacCormack predictor-corrector 
scheme, the first two operators Lj and L 2 are the two-step ex- 
plicit and dimensionally split method that MacCormack and 
Warming‘S proposed and Rizzi, et al.^ used for nonreacting 
flow. The last one, is an implicit modified Euler scheme 
that Lomax and Bailey‘S have analyzed and found to be 
second-order accurate. Note that this last operator accounts 
for chemistry only, while the first two treat only convection 
and no chemistry. The splitting procedure, therefore, neatly 
segregates these two rather different phenomena. 

Stability Analysis 

Without any approximation, Eqs. (1) may be expressed in 
nondivergence form as 


be 


be 




( 8 ) 


where c is the column matrix with elements e=(Cj,p, 
Wj,/?). Stability c ond itions can.be determined in the 
usual manner from an analysis of the eigenvalues of C/ and 
Cj. From L/ the condition on the step size 5/ ( =Ax“^) in the 
X- direction is 


5/ < min 

\ ^ 1 

(9) 


representing Eqs. (6) is stable if the necessary condition 

6<min(6y,a2) (11) 


is satisfied. 

Accuracy Analysis 

The sequence in Eqs. (6) is known to be accurate only to the 
first order of 6 because of the noncommutativity of the 
operators. A symmetric sequence, however, 

UlV=LjL2L,.L,L2L,Ul, (12) 

can be constructed by reversing the cycle of operators on each 
successive step, and in the next paragraph we will show that 
second-order accuracy is regained by so doing. 

Consider Eq. (8) where c now is the solution vector 
discretized over a computational mesh or net of individual 
points. The partial derivative terms C 2 (b/bx^) and Ci(b/ 
bx^) are represented by the matrices A and B and will later be 
approximated by the particular difference scheme that is 
chosen. So far, however, the full nonlinearity of Eq. (8) 
remains intact, but at this time A must be locally linearized by 
A= We-\-Q{b^) where IF is the Jacobian matrix [bA/bt]. 
Equation (8) may therefore be written as 

and further, if C/, C^, and W are constant overx^, the second 
derivative in x^ then becomes 

5^(5^-) = 1a + b.w,.. 

With the aid of this last assumption the accuracy of sequence 
in Eq. (12) will now be analyzed on Eq. (13), although with 
more labor it can also be proven that the symmetric sequence 
in Eq. (12) maintains the same accuracy as a second-order 
nonsplit method for the fully nonlinear problem as well. The 
explicit nonsplit MacCormack scheme when applied twice to 
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Eq. (13) yields the matrix equation 

«'>+' = |^/+26(A + B + W)+ (A^+B^+W^+AB 

+ BA + AW +WA + BW + WB)]e''+0(6^) (14) 

and is accurate to second order because it matches the Taylor 
series expansion of e through order 5^. In a similar way the 
sequence (10) has the matrix representation 

(/+^w)(/+5B + yB)(/+6A + y/l^)e'' 
or 

e'’+'=|^/+6(A + B + W + y (A2+5^+W^+2BA+2WB)je'' 

+ 0(6") (15) 

and because of the noncommutativity of the matrices the 
terms multiplying 6^ do not equal those in Eq. (14). Hence, 
the sequence is accurate only to the first order. Reversing the 
cycle for the next step, however, yields 

r 6^ 

c^+^=|^/+5(W + B + A)+y (W^ + B^+A^+2BW + 2AW 
+ 2AB)j[/+6(A + B + W)+y (/l^+B^+ B"^+2BA 

+ 2WBje"+ 0(6^) 

= |^/+26(A + B + W)+ -1-y- (A^+B^+W^+AB 

+ BA + AW + WA + BW + WB)je"+0(6^) (16) 

which is identical to the unsplit expression in Eq. (14) and 
*therefore the sequence (12) of mixed explicit and implicit 
operators in Eq. (6) is accurate to second order in space. 

Multiple Chemical Steps 

Splitting the difference scheme does not degrade its order of 
accuracy, instead some advantages do accrue. Not only does it 
simplify the programming but MacCormack and Warming*^ 
have pointed out that splitting the spatial dimensions in- 
creases the computational efficiency by allowing larger in- 
tegration steps. The greatest advantage, however, is gained 
from separating the fluid dynamical phenomena from those 
of chemical kinetics because each of these has intrinsically dif- 
ferent rates which they proceed. The convection process, for 
example, is slowly varying in space so that a relatively coarse 
step 6 can accurately resolve its variation in x\ whereas in 
comparison, the chemistry is much more rapidly changing and 
cannot be determined accurately through use of the same 
coarse step 5 even though the calculation is stable for any 6. 

To increase the resolution in this case we use a self- 
controlling increment routine that takes several fractional 
steps of 5,. for each coarse one of 5 and that, following Lomax 
and Bailey, checks at each step whether the last integration 
has produced a change in for the next partial step, and 
diminishes it if it has, so that in effect it always tries to ad- 
vance with a step size that yields a 10*7o change in the solution. 
From our experience in the cases presented here, the chemistry 
resolution typically requires five intermediate steps for each 
one for the transport process, thus severely increasing the 
computational time. Some workers have coped with this 
burden by solving the chemical rate equation [ first row of Eq, 
(4)] separately from the others. In their approach they ad- 



Flg. 3 Depiction of the variable-step procedure for L^. Several of the 
smaller chemical steps 6^. match one single larger conveclive slep 6. 

vance the gas dynamic equations at the usual CFL condition 
and then, with mixed explicit and implicit differences, return 
to integrate the rate equation using a variable step size similar 
to ours. 

While this intuitive procedure seems to give satisfactory 
results, it is difficult to determine its order of accuracy and 
stability criteria. Segregating the fluid and chemical 
phenomena by using the split operators in Eqs. (6) in the 
sequence 

E 

=L,(5)L,(5)YlLciK)L2mL,(5)Ul„ (17) 

( = r 

employing various chemical steps 5^ so that 

E 

Yi 6 . =26 

f = l 

is a more theoretically sound way of dealing with their dif- 
ferent rates. First it is more efficient because it allows all of 
the fluid to be convected by operators in Eqs. (6a) and (6b) in 
one single, large increment 6, and then remain idle during the 
time the chemical operator L, advances repeatedly at the 
smaller rale in order to “catch up” to the flow at 5 (see Fig. 
3). While in effect this is similar to what other methods do^’^ 
it is not identical because those methods advance the chemical 
rate equation which contains some convection terms that are 
differenced explicitly. Our L, contains no convection terms 
and hence eliminates the need for computing these differences 
as well as applying their boundary conditions at all in- 
termediate points. Furthermore, in contrast to the chemical 
rate equation, the split operator can be differenced easily 
with one implicit scheme that is unconditionally stable. Split- 
ting avoids the possible danger of restricted stablity that 
Lomax and Bailey‘s warn us can arise with use of mixed dif- 
ference expressions. 

Computational Mesh Construction 

In order for the sequence just described to advance the 
solution, the method must first construct the coordinate mesh 
network over which to integrate the equations. The coordinate 
system (x^ x*’), which may be nonorthogonal and 

unequispaced, is formed by three families of intersecting sur- 
faces (Fig. 4). The first one constructed is an arbitrarily 
oriented surface = constant at the n-\-I step, constrained 
only by the CFL condition on 5, and it intersects a second 
family = constant that at each extremity coincides with the 
shock** and the body surface and varies in some prescribed 
manner in between. The last family, x' = constant, is com- 
posed of simple planes and intersects the second at a number 
of given angular increments of <f>. This new marching surface 


**How to determine the location of the shock wave at the new step 
is explained in the next section. 
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MARCHING SURFACE 



Fig, 4 Definilion of general coordinates by the In- 

tersection of three surfaces, 

INTERMEDIATE MARCHING FINAL 



Fig. 5 Illustration of position, orientation, and vertex angle of 
conical data surfaces for successive steps. 


at /? + / can be completely general, and for the cases presented 
here we have chosen it to be circular conical because this 
requires at every step only one parameter for each of the three 
aspects of general motion. Translation of the cone is ac- 
counted for by the location x of its apex along the body axis, 
rotation by the angle between its axis and the body^s, and, 
lastly, dilatation by its vertex angle 0, as illustrated in Fig. 5. 
The second series of surfaces, of which the first and last are 
the body and shock, partition each cone into a number of con- 
tiguous annular segments. Each of these is then intersected by 
a plane rotated about the cone axis in specified angular in- 
crements, The resulting intersections delineate the coordinate 
and are straight lines or rays and complete the definition of 
the coordinates , x^, x\ As the solution proceeds down- 
stream, X advances while 0 increases monotonically toward 
90 and approaches zero so that at the final step the surface 
x^ - constant degenerates into a plane normal to the body 
axis. 

Shock Location and Boundary Conditions 

The intersections of the coordinates x^ and x^ lying in the 
marching surface x constant map out a mesh network of 
small quadrilateral ceils in such a way that the innermost row 
lies on the body and the outermost row is aligned with the 
shock. This method, termed mesh aligning, has the following 
features. For each cell at the outer edge of the mesh, the 
pressure Ps just downstream of the bow wave specifies the in- 
clination Q( the shoc k thro ugh the steady shock 

relations. The surface whose local slope matches this at each 
cell is constructed, and it becomes the outer edge of the mesh 
for the next step. Although the pressure p, at the shock is 
unknown it can be calculated (see Fig. 6) from the interior 



Fig. 6 Determination of the bow shock surface using a characteristic 
relation. 

values Pa by solving for the two unknowns, p^ and 
simultaneously and one local characteristic equation valid in 
the plane defined by the freestream and ij^ . The solution at 

each cell is carried out rapidly by an iterative procedure. 

Thus at the outermost row of cells becomes identical to 
the outward normal rj^ and to ^he body surface. This 
situation simplifies imposition of the conditions at the outer 
and inner edges of the overall mesh, referred to as entrance 
(shock wave) and streamline (body) boundaries. Along the en- 
trance boundary the flow variables are held fixed at their 
supersonic freestream values. Without modification, the dif- 
ference operators (6) then conserve mass and momentum 
across the discontinuity. This is precisely the condition used to 
derive the Rankine-Hugoniot relations, however, so that in ef- 
fect, application of the difference operators across the shock 
reduces numerically to the analytic jump conditions. Across 
cell faces that lie on the body, no transport is allowed. The 
only nonzero variable actually needed at that point in the 
calculations is the pressure which is related to the flow 
properties in adjacent interior cells by a characteristic relation 
in a manner similar to that at the shock wave. 

Results 

Computed results for reacting airflow about a general body 
traveling at 6.54 km/sec, Mach number M* =21.7, and 41° 
angle of attack are presented. The ambient freestream 
pressure is p^ - 106 dynes/cm, the temperature is T'* =236 K, 
and the body is a smooth three-dimensional configuration 
described by a series of third-degree polynomials in plan and 
profile views and ellipses in cross section. Between the body 
and shock the mesh holds nine cells and around the body in 
the meridional direction there are seventeen. 

The trace of the bow shock in the plane of symmetry, along 
with cross-sectional views at three axial locations, are 
displayed in Fig. 7, and clearly illustrate the large asymmetry 
of the flow at this high angle of attack. At cross section B-B 
the wing of the spacecraft is beginning to protrude out from 
the fuselage and grows steadily larger downstream. Figure 8 
presents the radial distribution of atomic oxygen Cq and nitric 
oxide Cno both the windward and leeward parts of the sym- 
metry plane at the three designated axial stations. The high 
temperature near the windward portion of the shock wave 
sharply enhances the concentration Cq of atomic oxygen 
which then quickly levels off nearer to the body. On the lee 
side, however, the shock temperature is much lower and the 
main rise in Cq occurs halfway between the shock and body, 
presumably due to convective transport from the windward 
region. For all axial positions, though, the same piaximum 
value is reached at the vehicle surface suggesting that the flow 
is probably frozen there. The concentration c^q of nitric 
oxide also increases dramatically to a maximum directly 
behind the windward portion of the bow wave and then falls 
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Fig. 7 Bow shock shape for reacting flow past a three-dimensional 
body (raveling at velocity #'qo= 6.54 km/sec, Mach number 
= 21.7, temperature =236 K, pressure = 106 dynes/cm and 
angle of attack = 4P . 

rapidly to a minimum at the spacecraft. On the lee side c^o is 
small at the shock, and the maximum is achieved roughly mid- 
way to the vehicle. The overall nature of the species 
distribution is typified in Fig. 9 by the circumferential 
variation of Cq and c^Jo along the shock and body in the plane 
Zax = 16.6 m. Since the bow shock is hottest on the windward 
side, the production of both O and NO is also largest there. 
Approaching leeward the shock becomes progressively 
weaker, and the concentrations of both O and NO fall 
precipitously. The same trend, however, does not occur at the 
body. Although Cno attains a maximum at about 70“ where 
the wing protrudes, Cq remains curiously constant at a 
maximum value around the entire vehicle even though its 
production at the lee region of the shock is very low. This in- 
dicates that local chemical production plays no role here. 
Rather, the convection process predominates; that is, atomic 
oxygen is produced at the windward part of the bow shock 
and is swept around to the lee side of the body. The vehicle is 
in an environment of greatly dissociated oxygen. 

For exactly the same flow as in Fig. 7, except that now a = 
30“ , two temperature plots in the windward symmetry plane 
are presented in Fig. 10 and compared with the results from a 
method of characteristics computation by Rakich et al.^ that 
are denoted by solid dots. The first part of Fig. 10 presents the 
axial variation along the body surface and the second at the 
radial distribution between body and shock at ^ov=15-7 m. 
Both indicate that temperature is a rather slowly changing 
property, usually in the range between 4000 K and 6000 K, 
although interestingly there is a sharp rise of about 1000 K 
near the body. Rakich et al. found a similar effect and at- 
tributed it to an “entropy layer” around the body. 

Concluding Remarks 

The method presented demonstrates that nonorthogonal 
curvilinear coordinates can be useful for computing practical 
fluid problems and in fact are the underlying basis of the in- 
tegral or finite-volume form of the governing equations. 
Unlike other marching methods, the finite-volume approach 
offers a wider choice of surfaces upon which the initial con- 
ditions may lie, and it allows the solution to advance in an ar- 





Fig. 8 Rudial distribution of (he chemical species O and NO in the 
windward and leeward sides of (he symmetry plane for the flow in Fig. 
7. a) Radial distribution of O; b) Radial distribution of NO. 


bitrary direction while at the same time it simplifies ap- 
plication of the body boundary conditions. Furthermore, this 
procedure lends itself to a new scheme of aligning the mesh 
with the bow shock, which is simple to implement, accurate, 
and consistent with the interior flow. Using the fractional step 
routine does not lower a difference scheme’s order of ac- 
curacy, but rather several benefits are gained. Split difference 
operators prove to be very appropriate for handling the di- 
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ferential distribution 
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Fig. 10 Temperature plots in the windward symmetry plane for flow 
identical to that in Fig, 7 except that now « = 30*; a) Axial variation 
along the body surface; b) Radial variation between body and shock. 


verse nature of the chemical kinetics and gas-dynamics 
phenomena. In particular, they provide a natural way of 
mixing explicit and implicit difference schemes that may use 


differing step sizes, and at the same time ease the com- 
putational burden. Results from this method are in accord 
with those from a characteristics procedure. Furthermore, a 
check on the conservation of total oxygen and nitrogen atoms, 
expressions independent of the calculations, reveals an overall 
error of only about 1%. The computation for the case presen- 
ted terminates at about 20 m because the shock wave 
generated by the leading edge of the wing collides with the 
bow wave and appears to produce a small pocket of subsonic 
flow at that point. The space-marching finite-volume method, 
of course, is restricted to supersonic flow and another 
procedure would be required to compute this phenomenon. 
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